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Abstract 

The hydrodynamic approach to a continuum mechanical description of granular behavior is re- 
viewed and elucidated. By considering energy and momentum conservation simultaneously, the 
general formalism of hydrodynamics provides a systematic method to derive the structure of con- 
stitutive relations, including all gradient terms needed for nonuniform systems. An important 
input to arrive at different relations (say, for Newtonian fluid, solid and granular medium) is the 
energy, especially the number and types of its variables. 

Starting from a careful examination of the physics underlying granular behavior, we identify the 
independent variables and suggest a simple and qualitatively appropriate expression for the gran- 
ular energy. The resultant hydrodynamic theory, especially the constitutive relation, is presented 
and given preliminary validation. 

PACS numbers: 
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I. INTRODUCTION 



When unperturbed, sand piles persist forever, demonstrating in plain sight granular me- 
dia's ability to sustain shear stresses - an ability that is frequently considered the defining 
property of solids. On the other hand, when tapped, the same pile quickly degrades, to form 
a layer (possibly a monolayer) of grains minimizing the gravitational energy This is typical 
of liquids. The microscopic reason for this dichotomy is clear: The grains are individually 
(and ever so slightly) deformed if buried in a pile, which is what sustains the shear stress. 
When tapped, the grains jiggle and shake, and briefly loose contact with one another. This 
is why they get rid of some of their deformation - which shows up, macroscopically, as a 
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gradual lost of the static shear stress and a continual flattening of the pile. 

When sand is being sheared at a constant rate, both solid and fluid behavior are operative. 
First, the grains are being deformed, increasing the shear stress as any solid would. Second, 



the same shear rate also provokes some jiggling, just as if the grains were lightly tapped. [43 1 
This leads to a fluid-like relaxation of the shear stress - the larger the shear rate, the stronger 
the jiggling, and the quicker the relaxation. Note the reason why loading and unloading give 
different responses (called incremental nonlinearity l|,|2|): When being loaded, the solid part 
of granular behavior increases the stress, while the fluid part decreases it. During unloading, 
both work in the same direction to reduce the stress. 

This entangled behavior, we suspect, lies at the heart of the difficulty modeling sand 
macroscopically. In addition, there is a "history-dependence" of granular behavior that, 
being experimentally obvious but conceptually confused and ill-defined, further perplexes 
the modeler. Obviously, if sand can be characterized, as do other systems, by a complete 
set of state variables, any history-dependence only indicates that the experiments were run 
at different values of these variables. This is what we believe happens. 

The hydrodynamic theory is a powerful approach to continuum-mechanical description 
(or macroscopic field theory), pioneered by Landau jsj and Khalatnikov [4] in the context of 
superfluid helium. Bei considering energy and momentum conservation simultaneously, and 
coining both with thennodyna,™ con S1 de r at 10 n S , this approach is capable of cogeatiy 
deducing, among others, the proper constitutive relation. Hydrodynamics [5| has since been 
successfully employed to account for many condensed systems, including liquid crystals 
0,0], superfluid 3 He js 10] , superconductors 11-13], macroscopic electro- magnetism jl4- 
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16| and ferrofluids 17H20L Transiently elastic media such as polymers are under active 



consideration at present 
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Two steps are involved in deriving the theory hydrodynamically, the first specifies the 
theory's structure: Being a function of the state variables, the energy itself is not indepen- 
dent. Nevertheless, the form of the energy density w(s,p) is left unspecified in this first 
step, and the differential equations are given in terms of the energy density w, its variables 
and conjugate variables. [Conjugate variables are the derivatives of the energy with respect 
to the variables, say temperature T(s,p) = dw/ds and chemical potential p(s,p) = dw/dp 
for s,p, the entropy and mass density]. In a continuum theory, a number of transport co- 
efficients [such as the viscosity i](s,p) or the heat diffusion coefficient n(s,p)] are needed 
to parameterize dissipation and entropy production. Neither is their functional dependence 
specified. 

A theory is unique and useful, of course, only when its energy and transport coefficients 
are made specific, in a second step. This division is sensible, because the first step is 
systematic, the second is not. The first starts with clearly spelt-out assumptions based on 
the basic physics of the system at hand, which is followed by a derivation that is algebraic 
in nature, and hence rather cogent. The second step is a fitting process - one looks for 
appropriate expressions, by trial and error, for a few scalar functions that, when embedded 
into the structure of the theory, will yield satisfactory agreement with the many experimental 
data. 

Starting from the physics of granular deformation and its depletion by jiggling, we have 
identified the variables and derived the structure of the equations governing their temporal 
evolution 24J |. calling it GSH, for granular solid hydrodynamics. But our second step is not 



yet complete, and some proposed functional dependencies are still tentative. The expression 
for the energy appears quite satisfactory, but our notion of the transport coefficients is still 
vague. Our final goal is a transparent theory with a healthy mathematical structure that 
is capable of modeling sand in its full width of behavior, from static stress distribution, via 



elastoplastic deformation [25|, [26(, to granular flow property at higher velocities [27H30]. 
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II. GRANULAR STATE VARIABLES 



In this section, we determine the complete set of granular variables starting from the 
elementary physics of granular deformation and its depletion by jiggling. 

A. The Elastic Strain 

If a granular medium is sheared, the grains jiggle, roll and slide, in addition to being 
deformed. Only the latter leads to a reversible energy storage. Therefore, the strain = 
Uij + Pij has two parts, the elastic and plastic one, with the first defined as the part that 
changes the energy. Hence the energy density w(uij) is a function of the elastic strain u^, 
which alone we identify as a state variable. For analogy, think of riding a bike on a snowy 
path, up a steep slope. The rotation of the wheel, containing slip and center-of-mass motion, 
corresponds to the total displacement d. The gravitational energy w(d t ) of the cyclist and 
his bike depends only on the center-of-mass movement d t , the "elastic" or energy-changing 
portion here. And the gravitational force on the center of mass is f g = —dw/ddt- Similarly, 
the elastic stress is nij = —dw(uij)/duij. When grains jiggle, granular deformation relax, 
hence 

d t Uij = Vij - Uij/r, (1) 

with the usual elastic term = |(VVu j + V/Ut), and a relaxation term —Uij/r that accounts 
for plasticity. [Note because the total strain obeys dtSij = vy, the evolution of the plastic 
strain p^ = —u^ is also fixed by Eq (pQ), and given as d t pij = u^/t] To understand how 
plasticity comes about, consider first the following scenario with r = constant. If a granular 
medium is deformed quickly enough by an external force, leaving little time for relaxation, 
J (uij/r) dt ~ 0, we have ity ~ = J Vijdt and = right after the deformation. The 
built-up in elastic energy and stress ir^ is maximal. If released at this point, the system snaps 
back toward its initial state, as prescribed by momentum conservation, d t (pvi) + Vj7Ty = 0, 
displaying an elastic, reversible behavior. But if the system is being held still {d t Eij = v»j = 
0) long enough, the elastic strain will relax, d t u^ = —Uij/r, while the plastic strain grows 
accordingly, d t Pij = u^/t. When vanishes, elastic energy w(uij) and stress 7Ty are also 
gone, implying dt (pVj) = 0. The system now stays where it is when released, and no longer 
returns to its original position. This is what we call plasticity. 
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However, 1/r is not a constant in sand: It grows with the jiggling of the grains (as the 
deformation is lost more quickly) and vanishes if they are at rest. If we quantify the jiggling 
by the associated kinetic energy, or (via the gas analogy) by a granular temperature T g , we 
could account for this by assuming 1/r ~ T g . 

As discussed above, a shear rate would jiggle the grains, giving rise to T g . For a constant 
rate, an expression of the form T g ~ y/VifvTj = || v s|| is appropriate [see Eq (Q3J) below]. 
Inserting 1/r = A||v s || (with A the proportionality coefficient) into Eq ([I]), we obtain the 
rate-independent expression, d t Uij = vy — Atty||v s ||. Being a function of ity, the stress 
TTij(uij) therefore obeys the evolution equation, 

dt-Kht = MkUjdtUij = Mkujiyij - Aiiy 1 1 v a 1 1), (2) 
M mj = diTki/duij = d 2 w/du ij du ki , 

which clearly possesses the structure of hypoplasticity , a state-of-the-art engineering 
model originally adopted because sand is incrementally nonlinear, and responds with differ- 
ent stress increases depending on whether the load is being increased (vy > 0, ||v s || > 0) 
or decreased (vy < 0, ||v s || > 0). It is reassuring to see that the realism of hypoplasticity is 
based on the elementary physics that granular deformation is depleted if the grains jiggle; 
and it is satisfying to realize that the complexity of plastic flows derives from the simplicity 
of stress relaxation. 

Under cyclic loading of small amplitudes, because the shear rate is not constant, T g 
oscillates and never has time to grow to its stationary value of T g ~ ||v s ||. Therefore, the 
plastic term u^/r ~ T g remains small, and the system's behavior is rather more elastic than 
rendered by Eq (j2J). 



The complete equation for ity is in fact somewhat more complex, 44] 



d t Uij = (1 - a)vij - m*-/t - u<u Sij/n, (3) 
1/r = XT g , 1/n = XiT g , (4) 

where u*j is the deviatoric (or traceless) part of and d t = d% + v^Vfe. The modifications 
are: (1) The relaxation time for u*j and uu are different. (2) A shear rate Vy yields an 
elastic deformation rate d^tty that is smaller by the factor of (1 — a). 

In contrast to strain relaxation ~ u^/t that is irreversible, a accounts for reversible pro- 
cesses (such as rolling). Without relaxation, elastic and total strain are always proportional, 



and for say a = 2/3, Uij is a third of £y. Circumstances are then reversible and quite analo- 
gous to a solid - aside from the fact that one needs to move three times as far to achieve the 
same deformation. So the physics accounted for by a is akin to that of a lever. [This is also 
the reason why the stress, or counter- force, is smaller by the same factor, see Eq (|19p .] Note 
since any granular plastic motion such as rolling and slipping, be it reversible or irreversible, 
become successively improbable when the grains are less and less agitated, we expect 

a(T g ) -» 0, for T g -» 0, (5) 

implying granular media are fully elastic at vanishing granular temperature. 



B. Mass, Entropy and Granular Entropy 

The energy density Wq(s, p) of a quiescent Newtonian fluid depend on the entropy density 
s and mass density p, both per unit volume. Defining the temperature and chemical potential 
as T = dwo/ds\ p and p = dwo/dp\ s , we note that they can be computed only if the functional 
dependence of wq(s,p) is given. The pressure, a prominent quantity in fluid mechanics, is 
also a conjugate variable, as it is given by P = dw/dv at constant sv, where v = 1/p is 
the specific volume, w = w v the energy per unit mass. Again, P is given once w (s,p) is. 
(Note it is not independent from p and T, since it may be written as P = — u>o + Ts + pp.) 

The conserved energy w depends also on the momentum density Qi = pvj, and is generally 
given as w = wo + g 2 /2p. So the complete set of variables is given as s,p and gt, and the 
hydrodynamic theory of Newtonian fluids consists of five evolution equations for them. Being 
a structure of an actual theory, these equations contain w , P, also T, p, Vj = dw/dgi. They 



are closed only when wo is specified. 45] 



In continuum-mechanical theories, the entropy s is not always given the attention it 
deserves. The basic facts underpinning its importance are: The conserved energy w is, in 
equilibrium, equally distributed among all degrees of freedom, macroscopic ones such as p, g^ 
and microscopic ones such as electronic excitations or phonons (ie, short wave length sound 
waves) . The entropy s is the macroscopic degree of freedom that subsumes all microscopic 
ones (typically of order 10 23 ), and accounts for the energy contained in them. Off equilibrium, 
energy is more concentrated in a few degrees of freedom, typically the macroscopic ones. 
The one-way, irreversible transfer of energy from the macroscopic to the microscopic ones 
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- in fluid mechanics from p, gi to s - is what we call dissipation, and the basic cause for 
irreversibility. A proper account of dissipation must consider the variable s, its conjugate 
variable T, and the entropy production R [with R/T denoting the rate at which entropy 
is being increased, see Eq (jUJ)]. This remains so for systems (such as granular media) that 
typically execute isothermal changes. 

The energy density of a solid depends on an additional tensor variable, the elastic strain 
Uij = Eij, which in crystals is very close to the total strain. The associated conjugate variable 
Ttij = —dwo/duij is the elastic stress - where linear elasticity, or -k^ ~ u^, represents the 
simplest case. The hydrodynamic theory of solids consists of eleven evolution equations, for 
the variables s, p, g^, u^, which in their structure contain the conjugate variables T, p, Vj, 7Ty. 

Displaying solid and liquid behavior, granular media have the same variables - in addition 
to the one that quantifies granular jiggling, for which a scalar should suffice if the motion 
is sufficiently random. We call it granular entropy s g , and define it to contain all inter- 
granular degrees of freedom: the stochastic motion of the grains (in deviation from the 
smooth, macroscopic velocity) and the elastic deformation resulting from collisions. We 
divide all microscopic degrees of freedom contained in s into the {46] inner- and inter-granular 
ones, s — s g and s g , with the conjugate variables T = dw /d(s — s g ) and T g = dw /ds g . 
Equilibrium is established, when both temperatures are equal, and s g vanishes. (There 
are overwhelmingly more inner than inter granular degrees of freedom. When all degrees 
have the same amount of energy, there is practically no energy left in s g .) The equilibrium 
conditions are: 

s g = 0, T g = T g -T = 0. (6) 
As zero is the value s g invariably returns to if unperturbed, it is an ener gy minimum. 



Expanding the s 9 -dependent part of the energy u? 2 = w — w(s g = 0), we take [471 ] 



w 2 (s,p, s g ) = s g /(2pb), T g = dw 2 /ds g \ s = Sg/pb, (7) 

with b(s,p) > 0. So the twelve independent variables are: s, s g , p, gi,Uij and the hydrody- 
namic theory consists of evolution equations for them all, of which six are given by Eq ([3]). 
The rest will be given in section IHIi These equations will contain wq and the conjugate 
variables: T, T g , p, v*, 7Ty, also the pressure, given as 

P T = —dwo/dv = -w + pp + sT + Sgfg, (8) 
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with the derivative taken at constant sv, s g v and Uij. As we shall see in Eq (I24p . this is the 
pressure that accounts for the contribution of agitated grains. 



C. History Dependence and Fabric Anisotropy 

Finally, some remarks about the special role of the density in granular behavior. First, 
it is quite independent of the compression ug£. Plastic motion rearranges the packaging 
and change the density by up to 20%, without any elastic compression. Second, the local 
density only changes if there is some jiggling and agitation of the grains, T g ^ 0. Even 
when non- uniform, a given density remains forever if the grains are at rest. So, if a pouring 
procedure produces a density inhomogeneity, this will persist as long as the system is left 
unperturbed, providing an explanation for the history dependence of static stress distribu- 
tion. Sometimes, these density inhomogeneities have a preferred direction, say, a density 
gradient along x. With density- dependent elastic coefficients, the system will then mimic 
fabric anisotropy, displaying a stress-distribution reminiscent of an anisotropic medium - 
even when it consists of essentially round grains and the applied stress is isotropic. Our 
working hypothesis, given a preliminary validation in section HV A 3\ is that both effects are 
covered by density inhomogeneities. The static stress of a sand pile is calculated there and 
compared to experiments for two densities, the first uniform and the second with a reduced 
core density, which we argue is a result of different pouring procedures, being rain-like and 
funnel-fed, respectively. 



III. GRANULAR SOLID HYDRODYNAMICS (GSH) 



This section presents the remaining six evolution equations. They will be explained but 



not derived, see 24| for more details and the complete derivation. 
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A. Entropy Production 



The evolution equation for the entropy density s is 



8 t s + Vi(svi - nViT) = R/T, 
+ 7 f s 2 + /3(vr*.) 2 + 



(9) 
(10) 



Eq Q is the balance equation for the entropy s. It is (with R unspecified) quite generally 
valid, certainly so for Newtonian fluids and solids. The term svi is the convective one that 
accounts for the transport of entropy with the local velocity, and kW{T is the diffusive term 
that becomes operative in the presence of a temperature gradient. R/T > is the source 
term. It vanishes in equilibrium, and is positive-definite off it, to account for the fact that the 
conserved energy w always goes from the macroscopic degrees of freedom to the microscopic 
ones, w — > s. 

The functional dependence of R changes with the system. In liquids, R is fed by shear 
and compressional flows, and by temperature gradients jsj] , as depicted by the first line of 
Eq ( ITOj) . In equilibrium, we have Vy, VjT = 0; off it, the quadratic form with positive shear 
and compressional viscosity, rj, C > and heat diffusion coefficient, k > 0, ensures that the 
entropy s can only increase. In fact, the terms of the first line are, in an expansion of R, 
the lowest order positive ones that are compatible with isotropy. 

The second line of Eq (II Op . with 7, (3, Pi > 0, displays the additional dissipative mech- 
anisms relevant for granular media. As discussed in the introduction, a finite T g or n^, 
indicating some jiggling or deformation of the grains, will both relax and give rise to en- 
tropy production. Since granular stress 7Ty will not dissipate for T g = 0, we require /3, /3i — >- 
for T g 0. 

Being part of the total entropy, the granular entropy s g obeys a rather similar equation, 
though it needs to account for a two-step irreversibility, w — > s g — > s, the fact that the 
energy goes from the macroscopic degrees of freedom to the mesoscopic, inter granular ones 
of s g , and from there to the microscopic, inner granular ones of s, never backwards, 





.2 
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(11) 

(12) 
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Eq (|TT|) has the exact same form as Eq (J9j) , so do the first three terms of R g . But R g also has 
a negative contribution. The three positive ones, with rj g , ( g , n g > 0, account for w — > s g , 
how shear and compressional flows, and gradients in the granular temperature produce s g , 
the jiggling of the grains. The negative term —jTg accounts for s g — > s, how the jiggling 
turns into heat. There is the same term, though with negative sign, in R, because the same 
amount of energy arriving at s must have left s g . As emphasized, all transport coefficients 
rj, Tj g , C, ( g , K i K gi 7) Pi Pi are functions of the state variables (which may alternatively be taken 
as T, f g , p, iTei and ttJ = 7r*-7r*-). 

In the stationary and uniform limit, for R g = and VjT 9 = 0, macroscopic flows produce 
the same amount of granular entropy as is leaving, implying 

lT g = %v*v*- + C g v 2 u . (13) 

This is the relation employed to arrive at Eq 02), showing that hypoplasticity holds in the 
limit of stationary shear rates. Given a shear rate, part of its energy will turn into s g , 
which in turn will leak over to s. At the same time, some of the flow's energy will heat up 
the system directly, with the ratio of the two dissipative channels parameterized by rj/rjg 
and C/C g - In dry sand, r], ( are probably negligible und shall be neglected below - though 
they should be quite a bit larger in sand saturated with water: A macroscopic shear flow 
of water implies much stronger microscopic ones in the fluid layers between the grains, and 
the dissipated energy contributes to R. 

Finally, we consider the T 9 -dependence of r] g ,(g,1- Expanding them, 

V = Vo + ViTg, ( g = Co + ClTgi 7 = 7o + JiTg, (14) 
we shall assume r] , ( Q = 0, because 

• R g then stays well defined for T g — >■ 0, see Eq (TT2l) ; 

• Viscosities typically vanish with temperature; 

• This fits the Bagnold scaling; 

• For 70 ^> 7iT 9 and 70 JiT g , respectively, we have from Eq ( !T3|) . for vu = 0, 

f g = (Vi/lo) |v*-| 2 , T g = |v*-|. (15) 

This ensures the existence of an elastic regime at vanishing T g , see section HVB2I 
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B. Conservation Laws 



The three evolution equations left to be specified are conservation laws, for mass, energy 
and momentum, 

dtp + Vi(pvi) = 0, d t w + WiQi = -pvJ7i<f>, (16) 
dt(pVi) + ^i(cxij + pViVj) = -pVi0, (17) 

where <fi is the gravitational potential (on the earth surface, we have — Vj0 = Gi, the 
gravitational constant pointing downwards). Without specifying the fluxes Qi,(Tij, these 
equations are always valid, quite independent of the system, and express the simple fact 
that being locally conserved quantities (in the absence of gravitation), energy, momentum 
and mass obey continuity equations. The basic idea of the hydrodynamic theory is to 
require the structure of the fluxes Qi, to be such that, with the temporal derivatives of 
the variables given by Eqs (151191 ll|ll 6 |fTTI) . the thermodynamic relation 

d t w(s,s g ,p,gi,Uij) = (dw/ds)d t s + (dw/ds g )d t s g + (dw/dp)d t p 

+(dw/dg i )d t g i + (dw/du^dtUij 
= Td t s + fgd t Sg + p,d t p + Vid t gi - Kijd t Uij 

is identically satisfied, irrespective of io's functional form. This is a rather confining bit of 
information, enough to uniquely fix the two fluxes as 

Qi = (W + P T )Vi + O-ijVj - KTViT - KgfgViTg, (18) 
C7y = (1 - a)7Ty + (P T - (gVujSij ~ TJgV^, (19) 

with P T given by Eq (ISL and v*- being the deviatory (or traceless) part of Vj.,-. (For de- 



tails of derivation see 



24|.) Although now specified to fit granular physics as codified in 
Eqs (!3f9][TT|) . these are still fairly general results, valid irrespective what concrete form w 
assumes. Moreover, they also nicely demonstrate the dependence on the number and types 
of variables: Eliminating s g , or equivalently, taking T g = 0, in Eqs (I8ll8|ll9p . one obtains 
the solid hydrodynamics. [48j Further eliminating by taking 7r.y = leads to the fluid 
hydrodynamics. 

Focusing on the plastic motion, the standard approach (especially the thermodynamic 
consideration by Houlsby and coworkers, [4oj]) employs the plastic strain = — Uij as 

11 



the independent variable. Although this starts from the same insight about plastic motion, 
the connection between elastic strain, stress and energy, so similar in solids and granular 
media, with formulas that hold for both systems, is lost - or at least too well hidden to be 
useful, see also the discussion in section II V C 1 1 

Enforcing a velocity gradient Vy, the rate of work being received by the system is 
-o-jjVjj = —[(l — a^ij+PTSi^Vij + lCgVitvu+rjgVijV^], see Eq ([18]). Of these, the terms in the 
first square brackets, being proportional to the velocity and hence odd under time inversion, 
are reactive; while those ~ v 2 in the second bracket are even and dissipative. Work received 
via an odd term will leave if its sign is changed by inverting time's direction; work received 
via an even term stays, as happens only with dissipative processes. The reappearance of the 
same factor (1 — a) as in Eq is not an accident, but required by energy conservation. 
If the same velocity leads to an elastic deformation that is smaller by (1 — a), then just as 
with a lever, the force counteracting this deformation = (1 — aju-ij + • • • is smaller by 
the same factor. 

This concludes the derivation of the structure of GSH, or granular solid hydrodynamics, 
given by Eqs (jgjgjl . (MT0]fTT]fT2l) and fll6fl7)18fl9p . 

IV. VALIDATION OF GSH 

The advantage of GSH is two-fold, its clear connection to the elementary granular physics 
as spelt out in the introduction, and more importantly, the stringency of its structure. It 
cannot be changed at will to fit experiments, without running into difficulties with gen- 
eral principles. The only remaining liberty is the choice of the functional dependence for the 
energy and some transport coefficients. As this implies much less wiggle room than with typ- 
ical continuum-mechanical models, any agreement with experimental data is less designed, 
"hand-crafted," and more convincing, especially with respect to the starting physics. 

In what follows, we shall fist examine granular statics, for a medium at rest, T g = 0, then 
go on to granular dynamics, with enforced flows or stress changes, and some accompanying 
jiggling, T g 7^ 0. An expression for the conserved energy w will be proposed that, in spite 
of its relative simplicity, reproduces many important granular features when embedded into 
GSH. As discussed above Eq ([H]), we divide w into three parts: the micro-, macro- and 
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mesoscopic ones, 



w 



w (s,p) + [w x {uij,pigi) +g 2 /2p] + w 2 (s g ,p). (20) 



The first 



49] accounts for the inner-granular degrees of freedom, all subsumed as heat into 
the true entropy s. We take wo = (E{s)/m)p, where E(s) is the energy of a grain, m 
its mass, and () denotes the average. The second consists of the contributions from the 
macroscopic variables of momentum density gi and the elastic strain u^, where w\ is given 
by Eq (|2"T|) below. The third, w 2 (s g ,p) of Eq (j7J), is further specified in section IIVB II 
It accounts for the inter-granular degrees of freedom, the mesoscaled, strongly fluctuating 
elastic and kinetic contributions. 



A. Granular Statics, T g = T 

Given an energy w%(uij), we can use the stress 7Ty(iiy-) = —dw%/duij and = \{ViUj + 
V jUi) to close the stress balance V jit^Ti) = pGi, and determine 71^- (r») with appropriate 
boundary conditions. As this is done without any knowledge of the plastic strain, we may 



3l| 



with some justification call this granular elasticity 

The relation Uij = |(VjZ7j + Vj£/j) remains valid because of the following reasons: In 
an elastic medium, the stressed state is characterized by a displacement field from a unique 
reference state, in which the elastic energy vanishes. Because there is no plastic deformation 
Uf, the total displacement is equal to the elastic one. Circumstances appear at first quite 
different in granular media. Starting from a reference state, a stressed one is produced 
by the displacement U + Uf, with typically Uf ^> C/j. Due to sliding and rolling, Uf 
is highly discontinuous, but U remains slowly varying, because the cost in elastic energy 
would otherwise be prohibitive. Fortunately, Uf is quite irrelevant: We have innumerable 
reference states, all with vanishing elastic energy and connected to one another by purely 
plastic deformations. As a result, we can, for any given displacement U + Uf, switch to the 
reference state that is separated from the original one by Uf, and to the stressed one by U. 
Now, the circumstances are completely analogous to that of an elastic medium. 



13 



1. Yield Surfaces 



An important aspect of granular behavior, in the space spanned by the variables, is the 
existence of yield surfaces. We take them to be the divide between two regions, one in which 
stable elastic solutions are possible, the other in which they are not - so the system must 
flow and cannot come to rest. A natural and efficient way to account for yield is to code it 



3l| - minimal if 



into the energy, a scalar. Given the stress balance, the energy is extremal 
convex and maximal if concave. Having the energy being convex within the yield surface, 
and concave beyond it, any elastic solution that is stable within the surface, will be eager to 
get rid of the excess energy and become unstable against infinitesimal perturbations beyond 
it. 



2. The Elastic Energy wi 



Our present choice for the elastic energy is 24|, I33M35] 



w 1 {p,u ij ) = BVA(2A 2 /5 + u 2 s /Z), (21) 

where A = —uu, u 2 = u*jU*j. The energy w\ is convex only for u s j A < \/2£, or equivalently 
its/ Pa < a/2/£ (where = |7r«, it 2 = tt^tt^), which coincides with the Drucker-Prager 



condition. 



50j Taking £ = 5/3 gives a friction angle of about 28°. We further take B 



Bq £>i(p) C(p,Uij), where Bq is a constant, and 



Bi = [(P-Pl p )/(P CP -P)r\ (22) 
2C = l + tanh[(A - A)/Ax]. (23) 

The coefficient B\ diverges for the "random closed-pack" density, p cp , and is convex only 
between p cp and the "random loose pack" density pi p . [p* t is a constant chosen to yield the 
right value for p ip with the relation p ip = (llp cp + 9p} p )/20.} It accounts for (1) the lack of 
elastic solutions for p < pe p , when the grains loose contact with one another; (2) the stiffening 
of granular elasticity with growing density, until it (as an approximation for becoming very 
large) diverges at p cp . 

With A , ki, k 2 , ks being constants, and A = kip — k 2 u 2 — k 3 , we have C = 1 for A <C A , 
and C = for A 3> Aq. It changes from 1 to in a neighborhood of Ai around Aq, destroying 
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FIG. 1: Yield surfaces as coded in the energy of Eqs (|21|22|23p . a function of the pressure, shear 
stress, and void ratio, (a): The virgin consolidation line, (b): The bending of the Coulomb yield 
line, as a function of e. (c): Combination of (a) and (b). 

the energy's convexity there. Taking A to grow with the density and fall with u 2 s limits 
the region of stable elastic solutions to sufficiently small A-values, reproducing the virgin 
consolidation curve and the so-called caps at varying void ratios e, see Fig [TJ 



3. Stress Distribution for Silos, Sand Piles and Point Loads 

Three classic silo, a sand pile and a granular sheet under a point load, are 

solved employing the stress expression derived from the energy of Eq (1211) . producing rather 
satisfactory agreement with experiments. 

Silos For tall silos, the classic approach is given by Janssen, who starts from the as- 
sumption that the ratio between the horizontal and vertical stress is constant, kj = a rr ja zz . 
Assuming in addition that a zz only depends on z, not on r, Janssen finds the vertical stress 
o zz saturating exponentially with height - a result well verified by observation. (He leaves 
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a rz and all three radial components: age, cr r g and a z e undetermined.) Having calculated <r zz , 
one needs the value of kj to obtain a rr , usually provided by kj ~ 1 — simp, with (p the 
friction angle measured in triaxial tests. This makes tp the only bulk material parameter 
in silo stress distributions. We shall refer to this as the Jaky formula, although it is also 
attributed to Kezdi. Being important for the structural stability of silos, this formula is 
(with a safety factor of 1.2) part of the construction industry standard, see eg. DIN 1055-6, 
1987. We believe this formula goes well beyond its practical relevance, that it is a key to 
understanding granular stresses, because it demonstrates the intimate connection between 
stress distribution and yield, a connection that has not gained the wide attention it deserves. 
Starting from Eq ( 121 p . we calculated 32) all six components of the stress tensor, verifying 
the Janssen assumptions to within 1%, and found the Janssen constant kj well rendered by 
the Jaky formula. 

Point Loads The stress distribution at the bottom of a granular layer exposed to a 
point force at its top is calculated {32 1 employing Eq ( 12TT) . without any fit parameter. Both 
vertical and oblique point forces were considered, and the results agree well with simulations 
and experiments using rain-like preparation. In addition, the stress distribution of a sheared 
granular layer exposed to the same point force is calculated and again found in agreement 
with experimental data, see 32( for more details and references. 

Sand Piles The fact that the pressure distribution below sand piles and wedges, 
instead of always displaying a single central peak, may sometimes show a dip, has intrigued 
and fascinated many physicists, prodding them to think more carefully and deeply about 
sand. Recent experimental investigations established the following connection: A single 
peak results when the pile is formed by rain-like pouring from a fixed height; the dip appears 
when the pile is formed by funneling the grains onto the peak, from a shifting funnel always 
hovering slightly above the peak. Employing Eq (l2Tj) to consider the stress distribution in 
sand wedges, we found the pressure at the bottom of the pile to show a single central peak if 
a uniform density is assumed. The peak turns into a pressure dip, if density inhomogeneity, 
with the center being less compact, is assumed. The two calculated pressure distributions 
are remarkably similar to the measured ones, see 31] . The nonuniform density, we believe, 
is a consequence of pile formation using the hovering funnel: Since the funnel is always just 
above the peak, the grains are placed there with very little kinetic energy, resulting in a 
center region below the peak that has a low density. Those grains that do not find a stable 
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position roll down the slope and gather kinetic energy. When they crash to a stop at the 
flanks, they compact the surrounding, achieving a much higher density. 



B. Granular Dynamics, T g 7^ T 

If a granular medium is exposed either to stress changes, or a moving boundary, the 
grains will flow, displaying both a smooth, macroscopic velocity, v, 7^ 0, and some stochastic 
jiggling, s g ~ T g 7^ 0. Then the following effects will come into play: First, the energy 
is extended by a s 9 -dependent contribution, w 2 (s s ,p), see Eq (JTj). Second, the transport 
coefficients of Eq (fT4"l) become finite. Most importantly, third, the relaxation times r, T\ of 
Eq 02]) are no longer infinite, implying the presence of plastic flows. 



1. The Sg- Dependent Part of the Energy 

Specifying the expansion coefficient b(p) of Eq (j7|) as b = b (l — p/p cp ) a , we find 

P T = apb,f 2 g {\ - plp cv ) a -\pl<lp cv ) (24) 

by employing Eq (jSJ). The density dependence of the expansion coefficient b(p) is chosen 
such that it reproduces the observed volume-dilating pressure contribution Pt ~ fi/ (Pc P — p) 



from agitated grains 
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However, we cannot take a = as it would imply a diverging 
granular entropy s g for p — > p cp . Therefore, we take a to be positiv but small, where a m 0.1 
appears appropriate. (Note that with w /p independent of p and w\/ p ~ A 2 5 - where A 
rarely exceeds 10 -4 - the respective density derivative and pressure contribution is zero and 
negligibly small.) 



2. The Hypoplastic Regime 

We may choose our parameters such that T g is small at typical velocities of elasto-plastic 
deformations, though large enough to cover both limits of Eq f)15p . Then the first term of 
Eq flTTTj) dominates, because all other terms (~ PT,rjg,( g ) are of order T 2 . Then we have 
dt&ij = (1 — ct)d t TCij = (1 — a)Mij k £d t Uij, with d t Uij given by Eq (J3]). Stress relaxation, the 
culprit producing irreversible plasticity, is a term ~ T g . For very slow shear flows and T g ~ 
I |vjj| | 2 [first of Eq f|T5|) ]. it is quadratically small and negligible. This is the elastic regime. At 
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somewhat faster shear flows, the relation T g ~ ||vf/|| [second of Eq ([175]) ] renders d t (Jij rate- 
independent, giving it the basic structure of hypoplasticity, Eq Q. Comparing this results 
to a state-of-the-art hypoplastic model, we found impressively quantitative agreement, see 
Fig[2J This is remarkable, because the anisotropy of these figures, determined essentially by 
Mijki, is a calculated quantity: M^t = d 2 Wi/duijduke, with w\ given by Eq (T2"TT) . 



3. The Butterfly Cycle 

Our last example for validation is not a direct comparison of ghd to some experimental 
data, but rather an examination of what ghd does, unforced and uncrafted, under typical 
elasto-plastic deformations. It is solved numerically for stress paths in the triaxial geometry 
(ie. a xx = <jyy, dij = for i ^ j, similarly for mj), including all energy terms given above, 
except C of Eq (1231) that is set to 1 (assuming the yield surface is sufficiently far away). All 
transport coefficients depend on T g as specified, but are otherwise constant, independent of 
stress and density. Also, all variables are taken to be spatially uniform, reducing a set of 
partial differential equations to ordinary ones in time. In spite of these major simplifications, 




FIG. 2: The change in strain d7 = (vn — V33)dt,de = — (2vn +V33)di for given stress rate starting 
from different points in the stress space, spanned by a s ,P, as calculated employing (1) GSM, the 
present theory (taking 1 - a = 0.22, rjr\ = 0.09, Cg/Vg = 0-33, Xy/Vg/l = H4), and (2) hpm, a 



typical hypoplastic model, see 



391 ] for more figures and details. 
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FIG. 3: Upper row: radial stress o\ versus axial stress 0-3, rescaled by Bqk~ 3 / 2 . Middle row: 
radial strain e\ = J v xx dt versus axial strain £3 = J v Z zdt. Lower row: e — eo (with eo the initial 
void ratio) versus shear strain e q = J(v 22 — v xx )dt, rescaled by v\K. The stress loads are isobaric 
for (a) and quasi-isobaric for (b,c); the cyclic amplitude is small for (a,b) and large for (c). The 
associated strain loci and void ratio are: sawtooth- like for (a), coil-like for (b), butterfly-like (or 
double- looped) for (c). [The large-amplitude, isobaric plot is quite similar to (c).] 

the results as rendered in Fig 3 display such uncanny realism that it seems obvious GSH has 
captured some important elements of granular physics. We consider a test with the stress 
given as 

p = pav + pampl ^ ^ f ^ ? g = g ampl ^ ^ ff + ^ _ ^ 

Numerical solutions were computed for isobaric test with P am P l = (ie. P = constant) 
and quasi-isobaric test, with p am P l << p av (ie. P constant). The results are shown 
in Fig they are obtained using the dimensionless parameters: k = VC1I1/ pb = 18257, 
(7o/7i) 2 (Pc P b 2 n 3 / 2 /B b ) = 1.07 x lO" 6 , Ai/A = 0.09, uf/2 = Vl /3Ci = 1, X^fnUll = 114, 
a = 0. The initial conditions are: eo = 0.68085 (or p = 0.94p cp ), Vij, T g , d t T g , d t p, 
d t Uij = 0. The averaged pressure P av = is P av = 70B k~ 3 / 2 , and the amplitude 

g ampl = ^ _ ai ig buiK 3/2 q ampl/Q Bo = 1Q for &nd 1Q Q for ^J) The f requency Q f 
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P ampl , g^P 1 is f — 12(7o/&p), and the phase lag between them is ip = 58°. 



C. Competing Concepts and Misconceptions 

Finally, we revisit two previous approaches to come to terms with granular behavior, 



40] , and granular statistical mechanics by Ed- 



granular thermodynamics by Houlsby et al 
wards et al 4l| . We shall compare GSH to both assuming at most superficial familiarity with 
them. Also, we refute some misconceptions that have become unfortunately widespread, es- 
pecially the one about energy not being conserved in sand [sic]. These are at best a nuisance 
in exchanges with referees; and at worst actual obstacles in the progress of our coming to 
grips with granular modeling. 

1. Granular Thermodynamics 

Although considerable work and thoughts have gone into applying thermodynamics to 



granular media and plastic flow, especially from Houlsby and Collins [401 ]. its basic points 
are clear and easy to grasp. Taking the entropy production as 

R = ;:,,(), i>i j (26) 

(where Pij denotes, as before, the plastic strain), it is obvious that the usual linear Onsager 
force-flux relation, d t pij ~ iiij, hence R ~ 7r?-, does not give a rate-independent R. Therefore, 
Houlsby, Collins and coworkers consider instead 

R = a/ XijkedtPijdtPke = {XijkidtPijdtPu) / \JXijktdtPijdtPkt, 

a rate-independent expression. Equating it to Eq ( 1261) . with = —dF/dpij, and F being 
the free energy density, one then solves for the plastic strain pij with a given F. One example 
gives d t Pij ^ on a yield surface, characterized by some components of ir^ being constant, 
and d t Pij = off it. 

GSH starts with the same R, but possesses the additional variable T g , for which T g ~ ||v s || 
frequently holds, see Eq ( |T3l) . The linear Onsager force-flux relation 

dtPij = PiTij with ~ T g , (27) 
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therefore suffices to yield an rate-independent R ~ Tg^. Note Eq (|27|) leads directly to the 
relaxation term: Because d t Uij + d t pij = Vy, we have d t Uij — Vy = —fiit^ = —Uij/r, with 
1/r ~ T g . (The last equal sign holds because 7Tjj,/3,T are all functions of w^, with (3,r as 
yet unspecified.) 

Summarizing, without the variable T g , Houlsby and Collins needed to go beyond the 
well- verified and -substantiated procedure of linear Onsager force-flux relation to maintain 
rate-independence, obtaining a plastic flow that is confined to the yield surface. In GSH, 
rate-independence arises naturally within the confines of linear Onsager relation, producing 
a plastic flow that is as realistic as hypoplasticity, and finite also off the yield surface. 



2. Granular Statistical Mechanics 



Generally speaking, it is important to remember that of all microscopic degrees of free- 
dom, the inner-granular ones are many orders of magnitude more numerous than the inter- 
granular ones. It is the former that dominate the entropy and any entropic considerations. 
When revisiting granular statistical mechanics, especially the Edwards entropy, it is useful 
to keep this in mind. 

Taking the entropy S(E, V) as a function of the energy E and volume V, or dS = 
(l/T)dE + (P/T)dV, the authors of [4l| argue that a mechanically stable agglomerate of 
infinitely rigid grains at rest has, irrespective of its volume, vanishing energy, E = 0, dE = 0. 
The physics is clear: However we arrange these rigid grains that neither attract nor repel 
each other, the energy remains zero. Therefore, dS = (P/T)dV, or dV = (T/P)dS = XdS*. 
The entropy S is obtained by counting the number of possibilities to package grains for a 
given volume, and taking it to be e s . Because a stable agglomerate is stuck in one single 
configuration, some tapping or similar disturbances are needed to enable the system to 
explore the phase space. 

In GSH, the present theory, grains are neither infinitely rigid, nor always at rest, hence 
the energy contains both an elastic and a s 9 -dependent contribution. [5l| And the question 
is whether granular statistical mechanics is a legitimate limit of GSH. We are not sure, 
but a yes answer seems unlikely, as both are conceptually at odds in two points, the first 
more direct, the second quite fundamental: (1) Because of the Hertz-like contact between 
grains, very little material is being deformed at first, with the compressibility diverging at 
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vanishing compression. This is a geometric fact independent of how rigid the bulk material 
is. Infinite rigidity is therefore not a realistic limit for sand. (2) As emphasized, the number 
of possibilities to arrange grains for a given volume is vastly overwhelmed by the much 
more numerous configurations of the inner granular degrees of freedom, especially phonons. 
Maximal entropy S for given energy therefore realistically implies minimal macroscopic 
energy, such that a maximally possible amount of energy is in S (or heat), equally distributed 
among the inner-granular degrees of freedom. Maximal number of possibilities to package 
grains for a given volume is a very different criterion. 

3. Energy Conservation 

Stemming ultimately from a loose vocabulary, some alleged difficulties to model sand are 
based on fallacies that need to be refuted here. 

The essential difference between granular gas and ideal (atomic or molecular) gas is that 
the particles of the first undergo non-elastic, dissipative collisions. As a result, their kinetic 
energy is not conserved, and the velocity distribution typically lacks the time to arrive at the 
equilibrium Gaussian form. Quantifying the kinetic energy as a granular temperature T g , 
it is therefore hardly surprising that the fluctuation-dissipation theorem (fdt), formulated 
in terms of T g , is frequently violated. These are sound results, obtained from a healthy but 
truncated model that takes the grains as the basic microscopic entity with no heat content. 
However, some of the further conclusions are deduced forgetting this simplification, rendering 
them patently absurd. These, and their [refutation in italic], are listed below: 

• As the energy is not conserved in sand, neither thermodynamics nor the hydrodynamic 
method are valid. [Only the kinetic energy dissipates in granular media, not the total 
energy. The latter, including kinetic, elastic and heat contributions, remains conserved 
- as it is in any other system. And only the conservation of total energy is important 
for thermo- and hydrodynamics.] 

• fdt, along with other general principles either derived from it or in its conceptual 
vicinity (such as the Onsager reciprocity relation) are all violated. [There are two 
versions of FDT, only the one given in terms of T g is violated, not the one in terms 
of the true temperature T. The latter is a general principle and always valid. For 
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instance, the volume fluctuation is given as (AV 2 ) = T(d 2 F/dV 2 )~ 1 , with F the 
associated free energy, for a copper block, a single grain, and a collection of grains. If 
the grains in the collection are jiggling, there is an extra contribution ~ T 2 in F, see 
Eq FK), that considerably increases the value of (AV 2 ) . The Onsager relation remains 
valid because the true fdt holds.] 

• The Onsager relation is also violated because the microscopic dynamics, the collision 
of the grains, is dissipative and hence irreversible. [The true microscopic dynamics is 
that in terms of atoms and molecules, the building blocks of the grains. Their dynamics 
is, as in any other system, reversible.] 
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